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ABSTRACT 

We consider positivity preserving property of first and higher order finite volume schemes 
for one and two dimensional compressible Euler equations of gas dynamics. A genera] frame- 
work is established which shows the positivity of density and pressure whenever the under- 
lying one dimensional first order building block based on an exact or approximate Rie- 
mann solver and the reconstruction are both positivity preserving. Appropriate limitation 
to achieve high order positivity preserving reconstruction is described. 
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1 Introduction 


To solve a scalar conservation law 


u t + div f(u) = 0 (1.1) 

with possibly discontinuous solutions, one usually studies the total variation stability (in one 
space dimension) or L°° stability (in multi space dimensions) of the numerical schemes. The 
schemes can usually be generalized to systems via local characteristic decompositions and 
usually work equally well numerically. However, no stability property can be automatically 
carried over to the nonlinear system case. For example, most second order TVD schemes, or 
even some first order monotone schemes, when generalized to compressible Euler equations 
of gas dynamics, do not always preserve the positivity of density and pressure. This may 
cause problems in practical calculations when the solution is near vacuum, for example in 
the computation of blast waves. 

When considering the compressible Euler equations of gas dynamics, a natural stability 
condition for the numerical approximation is positivity preserving for density and pressure. 
Another possible a priori bound is the maximum principle for the specific entropy (Tad- 
mor [10]), which seems extremely difficult to preserve for second and higher order schemes 
(Khobalatte and Perthame [5]). 

The positivity of density and pressure is already interesting because it allows one to derive 
a rigorous CFL condition (limitation of the time step to use). Since we consider conservative 
schemes, it also allows one to obtain a priori L l bounds on the density, momentum and energy. 
In this sense, positivity of density and pressure is a weak stability condition. Of course for 
these nonlinear systems, such a weak stability is not enough to assert the convergence of the 
method and estimates on the derivatives are usually needed for that purpose. Note however 
that for a large class of 2 x 2 systems in one dimension, weak stability (say in L°°) is enough 
to prove the convergence of the method. 

In this paper we will provide a general framework and illustrate by several examples the 



way to impose positivity of density and pressure for finite volume schemes, for one and two 
space dimensions and for first and higher order accuracy. Of course the first ingredient is a 
positivity preserving exact or approximate Riemann solver, such as Godunov, Lax-Friedrichs, 
Boltzmann type (Perthame [ 6 ]), and Harten-Lax-van Leer [3]. For a first order scheme, this 
is enough to guarantee positivity of density and pressure, even for two dimensional problems 
set on an arbitrary triangulation. However, for higher order finite volume schemes, the 
reconstruction must also satisfy such a property. We will show that only the nodal values 
needed for constructing the flux along a cell edge must be positive, in order to obtain a 
positive scheme. The reconstructed function, which is a piecewise polynomial and can be 
obtained in ENO spirit, needs not be positive everywhere. 

This paper is organized as follows: we first prove a positivity result and apply it to a 
third order scheme in one space dimension in Section 2. Then, in Section 3, we consider 
first order positive schemes in two space dimensions for arbitrary triangulations. The fourth 
section is devoted to second order schemes in two space dimensions. In the Appendix we 
recall why Lax-Friedrichs scheme is positivity preserving. 

2 Positivity in One Space Dimension and a Positive 
Third Order Scheme 

In this section we consider the one dimensional compressible Euler equations for perfect gas: 

Ut + F(U) x = 0, t>0,xeR ( 2 . 1 ) 

with 

U = (p,pu,E), E=^pu 2 + pe (2.2) 

F(U) = (pu,pu 2 + p, (E + p)u) , p=( 7 -l)pe (2.3) 

where p is the density, u is the velocity, E is the total energy, p is the pressure, e is the 
internal energy, and 7 > 1 is a constant (7 = 1.4 for air). 
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Let us recall for later purpose that for the pressure law in (2.3), the speed of sound is 


given by c— ^ 7(7 — l)e and thus the three eigenvalues of the system are u,u±c. 

We give conditions on numerical schemes in order to satisfy the positivity property, by 
which we mean that the resulting value of U should satisfy p > 0 (p™ > 0 on the discrete 
level) and p > 0 (p? > 0 on the discrete level), if the initial condition satisfies those positivity 
conditions. 

We first present the finite volume schemes under consideration, then we give a general 
positivity theorem, and finally we show how a third order reconstruction can be modified in 
order to satisfy the assumption of the general theorem. 

A general finite volume scheme can be written as 







(2.4) 


where n > 0 refers to the time step At and j € Z to the uniform space discretization of 
the size Ax (to simplify the presentation), and A = Z7™ are approximations to the cell 

averages of U in the cell Cj — (Xj_i,x J+ i) at time level n, and Uj+i are high order 

approximations of the nodal values U n (xj + 1 ) within the cells C, and C,+i, respectively. These 
values can either be evolved as independent degrees of freedom such as in the discontinuous 
Galerkin finite element method (see, e.g., [ 1 ]), or be reconstructed from the cell averages 
IT-. Let us recall that the general ENO philosophy [4] allows us to reconstruct, from the cell 
averages Z7™> a piecewise polynomial function U n (x ) which is high order accurate (r-th order 
if U n (x ) is piecewise (r — l)-th order polynomials), and conserves the local mean: 

4- [ u n (x)dx = Tf] (2.5) 

Ax Jo, 

The nodal values needed in scheme (2.4) can then be set to 



since the function U n (x) is discontinuous at the cell interface x J+ i. 


(2.6) 
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What remains to be explained for scheme (2.4) is the flux function h(lJ, V). This is 
assumed to be an exact or, in most cases, an approximate Riemann solver. In particular, 
h(U, V ) is a Lipschitz continuous function of both arguments, and is consistent with the 
physical flux F(U) in (2.3): h(U,U ) = F(U). We will further make the following 

Assumption 1: produces a one dimensional first order scheme which satisfies the 

positivity property under a CFL condition: 

A II(M + C )lloo ^ a o (2-7) 

i.e., if p” > 0 and p" > 0 for all j € Z, then the solution f/” +1 of 

Uj*' = Uj - A [*(£/"+, , Uf) - h(Ul )] (2.8) 

also satisfies /9™ +1 > 0 and > 0 under the CFL condition (2.7). 

We will call such a /*(-, ■) positivity preserving under the CFL condition (2.7). Examples 
of positivity preserving (approximate) Riemann solvers include Godunov, Lax-Friedrichs (see 
the Appendix), Boltzmann type [6] and Harten-Lax-van Leer [3]. 

Our first result is 

Theorem 1. Assume /*(♦,*) satisfies Assumption 1. If the reconstructed nodal values (/* L 

3 + 2 

have positive density and pressure for all j £ Z, then the full scheme (2.4) is positivity 
preserving under the CFL condition 

A ll(kl + c )lloo (2.9) 

where 0 < a < 1 is sufficiently small such that 

u"j - °(V- . +(/;_,) := Vj (2.10) 

3^2 3 2 

has positive density and pressure for all j 6 Z. 

Remark 1. Of course a — 0 works. However, since L are close to 77™, we can expect 

J'F 2 3 

that a is close to i. 
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Remark 2. Here we restrict ourselves to first order Euler forward in time. TVD type high 
order Runge-Kutta time discretization (Shu and Osher [9]) will keep the validity of Theorem 
1, as it was shown in [5]. 


Proof of Theorem 1: Let us introduce three “first order” schemes within the cell Cf. 


Vt 


u: 


u 


J + 5 


Vj - A 


11 U 
J 2 


A 

a 

h 

A 


ft (</+,, (/-J 


( 2 - 11 ) 


The three of them are of the type (2.8), with possibly £ in the place of A. Therefore under 
the CFL condition (2.9), the values f7 + , U~ and V all have positive density and pressure. 
We set 

f7™ +1 = ot{Uf + U~) + Vj (2.12) 


Then the linear combination of the above three equalities in (2.11) with weights (a, 1, a) just 
gives the scheme (2.4), thanks to the definition of Vj in (2.10). 

Finally, (2.12) implies that, by concavity, ZT^* 1 has also positive density and pressure. 
The pressure is indeed a concave function of (p, pu, E). 


□ 


Theorem 1 tells us that, if we use a positivity preserving approximate Riemann solver, 

then we only need to worry about positive density and pressure in the nodal values i , 

either from the reconstruction or from direct evolution. This can usually be achieved by a 

further limitation upon U ± L , such that positivity of density and pressure is enforced and 

. 7+2 

accuracy is preserved in smooth regions. We now show, as an example, how this can be 
done for a third order reconstruction. Let us consider a typical cell ( — 8,8) where 8 — ^r, in 
which we are given three cell averages (p, pu, E) with positive density and pressure. Let us 
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build three quadratic polynomials 


x 

p { x ) = Po + P\X + P2Y 

X 1 

u(x) = u 0 + UiX + w 2 — (2-13) 

X* 

p(x) = po + p\x + Pi ~2 

obtained by (i) fitting the average of ( p,pu , By- + ^-) to (p,pu, E) and (ii) fitting the nodal 
2 

values ( p,pu , + ^y)(±6) with the values obtained by a third order reconstruction, say 

from the ENO procedure. These might be associated to a negative density and pressure. 
Then we need an additional limitation. It can always be performed easily as shown in 

Proposition 2: Given the three quadratic polynomials in (2.13), we can always perform 
a limitation of their coefficients so that (i) the means of the three conserved variables are 
preserved, and (ii) the nodal values have positive density and pressure. 

Remark: The precise limitation, and the deduction of the appropriate coefficients in (2.13) 
are given in the proof below. 

Of course, combining Theorem 1 and Proposition 2, we obtain a positive preserving 
scheme. A similar construction for second order schemes can be found in [6] and [5]. 

Proof of Proposition 2: First of all, let us fix some notations. Being given the cell averages 
(p,pu, E ), we also define u and p through 

W = pu, E=^- + —^— (2.14) 

2 7 — 1 

Notice that we are given p > 0 and p > 0. 

First step: The limitation on the density is rather easy. Since p = po+P 2 ^, we can modify 
p\ without altering the conservation of mass. To enforce positivity of p(±$) is equivalent to 
impose a simple limitation on pj: 



But the right-hand-side of (2.15) could be negative. Therefore we first impose, for instance, 

-y > po > ~y-> p~Po + P 2 y (2.16) 

4 4 o 

Indeed, this implies that 

\p *\~2 ~ 3 l^ 2 1 = ~ Po\< < Po 

and now (2.15) can be imposed too. 

Second step: Next the conservation of momentum implies the following relation which 
gives u 0 in terms of U],u 2 and p(x): 

P , 82 8 \ to 1 7\ 

pu 0 - pu- p\U\ — - u 2 (p 0 — + p2 7^) \ ZA< ) 

And the conservation of energy gives, after some easy calculations 

p=( 7 -l) -^-p{u 0 - uf + u\ ^p 0 y + /»2 + u x u 2 p\— + P 2 ^l— +Po + P2g- (2.18) 

As in the density case, we will impose, in order to get a positive pressure at the nodes, a 
simple limitation on p\ 

|pi |<5 < Po + P2~ (2.19) 

and thus we need some limitations in order to guarantee that the right-hand-side of (2.19) 


is positive. As before we first impose the limitations on po,iti,u 2 : 

PO < \v (2.20) 

«i(7-l)(4 + l/ ’ ll fo + ' >2 To) <fij (2 ' 21) 

•4(7 -0(1,4 + 4) <15 (2 ' 22) 

This defines a unique u 0 through (2.17) and a unique p 2 through (2.18). For these values we 
deduce from (2.18), using (2.20): 


P _ 2 

6 < P ~3 P ° 

1 / ^2 ^4 \ £4 £61 1 / 

= (7-1) -7iP( U 0-u ) 2 +uUp 0 — + P2-) + P1U1U2— + p^ul—\ + -[P0+P2 


, r 2 / P P . t P\ 2 / 8* , . . 
< (7-1) u i + ^ 2 jo + j + U2 v 2 56 + ^ 
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and from (2.21)-(2.22) we deduce that po + P 2 y > 0 and thus (2.19) can also be imposed. 

This concludes the description of a possible set of limitations, which proves Proposition 2. 
Notice that these limitations (2.15), (2.16), (2.19), (2.20), (2.21) and (2.22) just amount to 
avoid gradients or second order derivatives of the order ^ (or larger) and thus they preserve 
the accuracy in smooth regions. 

3 First Order Positive Schemes in Two Dimensions 


We now consider the equations 


U t + div F(U) = 0 t > 0, x € R 2 


(3.1) 


where 


1 


U = (p,pu,E), uZR 2 , E = -p\u\ 2 + pe 


(3.2) 

F(U) = (pu,pu <g> u + pI,(E + p)u), p={^-\)pe (3.3) 

We consider finite volume schemes set on a triangulation C. The control volumes will be 
the triangles K 6 C. For each triangle K we denote by (/^)i <«<3 the length of its three edges 
( e A')i<a <35 with outward unit normal (^-)i<a<3- Finally, K(ca) will denote the neighboring 
triangle along e % and | A" [ the area of the triangle K. Then, we consider the scheme 

At _ 2 

up' = K- - j^i E MCkh, (3.4) 

for some (approximate) Riemann solver h(U, V , v) in the direction v. We recall some classical 


examples in the Appendix. The basic properties of h are now 


h,(U,V, u) = — k(V, U, — v), ( conservativity ) (3.5) 

h(U,[J,v) = F(U) ■ v ( consistency ) (3-6) 

Next, we also impose a positivity condition for the one dimensional solver obtained by 
fixing v (see the Appendix for examples). This is the same as Assumption 1 in the previous 
section, which now means that the solution U, a four component vector, of 

U = U-\[h{V,U,v)-h(U,W,v)] (3.7) 


8 



has positive density and pressure as soon as U , V , W do. Now, h is a four component vector. 
The main result of this section is 

Theorem 3: Let /*(•,•,•) satisfy (3.5)-(3.6) and the one dimensional positivity property As- 
sumption 1. Then the scheme (3.4) satisfies the positivity property under the CFL condition 

Y + c K ) < ao|A"|, f or all K (z C (3.8) 

c* — 1 


Proof of Theorem 3: We define 


U a = 


/« 

l K 


EL, 1*1 


//” 

0\ V K 


At 


% l a K \hm (a) ,U n K ,u%.) - h{U' K ,U n K ,v a K)] (3.9) 


Since 


E W u k< Vk,*k) = f(Uk) ■ E "k'k = 0 


a=1 


a=l 


we have = X^ =1 U a . And we conclude as in the one dimensional case. 


4 Second Order Positive Schemes in Two Dimensions 


We extend the result of the previous section to second order schemes, i.e., in which edge 
averages of the flux are approximated with second order accuracy. It has to be noted that, 
as usual, this means that the resulting scheme formally only has a first order truncation 
error. In the same way the first order schemes used in Section 3 formally only have zeroth 
order truncation error but it is well known now that they are convergent (see for example the 
preprints of Szepessy entitled, “Measure valued solutions to conservation laws with boundary 
conditions”, of Champier, Gallouet and Herbin entitled “Convergence of an upstream finite 
volume scheme on a triangular mesh for nonlinear hyperbolic equations”, and of Coquel and 
LeFloch entitled “The finite volume method on general triangulations converges to general 
conservation laws”). 

As for the one dimensional case, we give a general result assuming a positive reconstruc- 
tion. 
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We now consider an approximation of the two dimensional Euler equation (3.1) under 
the form 


A t 3 

W' =U' K - tttt £ M 


(4.1) 


or= 1 


where we use the notations of Section 3. We still assume that the approximate Riemann 
solver h satisfies the one dimensional positivity property Assumption 1. The main difference 
with the situation in Section 3, is that we now use second order approximations ^x{a) 

of U n at the center of the edge e^. U^’ a is such an approximation in the triangle K , (/£’?) is 
in the triangle K(a) neighboring K along the edge e%. These values can be obtained using 
a slope reconstruction, or an interpolation together with the interpretation of functions as 
piecewise constant in subcells as in Perthame and Qiu [7], or being evolved as independent 
degrees of freedom in the discontinuous Galerkin finite element method [2], In any case, we 
can assume that these values are conservative, i.e. that there exists real numbers (w/*-)i< 0 < 3 , 
such that 

= 1*1. ^>0; = w (4.2) 

(3 = 1 

For instance, when U^ a is the value, in the middle of the edge e£-, of a linear function 
in K , the coefficients are just \K\ times the bary centric coordinates of the mass center 
of K , with respect to the middle of the edges. Another example can be found in [7]. 

In order to state our result, let us introduce some notations: Denoting by Cj*- the middle 
point of the edge e^, we define 

3 a rux 

H = (4.3) 


a 1*1 

then, the triangle K is naturally divided into three subtriangles (A" Q ) 1<a < 3 obtained by 
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joining Yk to the three vertices of K. 



Fig. 1 Sub-triangles decomposition of a triangle A 

The unit outward normals of the sub-triangle K a are denoted by )i<o,'<3 where 

and the length of the edges (e^-’ Q )i< a '< 3 , where e^ a = e^-, of the sub-triangle 
K a , are respectively (/^-’°' )i<<*'<3 with / = l k- 
We can now state our main result: 


Theorem 4: Let h satisfy the one dimensional positivity property Assumption 1, and assume 
(4.2). Then the scheme (4.1) satisfies the positivity property under the CFL condition 


£ /r'MK-l + c a K ) < a 0 u a K 

a'—l 

Proof of Theorem 4: Mimicking the proof of Theorem 3, let us define f/£ +1, “ by 


where 


Trn+l,oi a 7 rn,c* \ ^ h ( 1 TI n,a 

Uk — U K ' U) K — Z\t 2_^ h \ U K ) U K ) l K 

0=1 


u% a ’ a = U*'° a) , U^ a ’ 0 = U^ for ft 


Adding the equalities (4.5) for a = 1,2,3 we obtain (4.1) with 

vT' = E 

Indeed, this is just a consequence of (4.2) and, for a ^ ft, 

h(uz° J ’,uz A ,v°/) = -h(u"/- a ,u n /,v^) 


(4.4) 


(4.5) 


(4.6) 


(4.7) 


(4.8) 


__ //^1« 
l K — l K 


(4.9) 
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Now, applying Theorem 3 to (4.5), we obtain that the f/£ +1,a satisfy the positivity prop- 
erty. By concave combination f/£ +1 satisfies it too and the theorem is proved. 


□ 


Notice that (4.2) is not essential. It allows us to simplify greatly the statement of Theorem 
4. If it is not satisfied, we require a two dimensional extension of the assumption (2.10) for 
the one dimensional case which could be quite difficult to check here. Finally, we refer to 
[7] for an example of a reconstruction which satisfies the positivity property, together with 
(4.2). There the control volume is a dual mesh and the positivity is proved using Boltzmann 
solvers. Our approach here could be extended to a dual mesh and the full scheme, using a 
Lax-Friedrichs or Godunov solver, satisfies also the positivity property. 

5 Conclusion 

We have considered how to preserve the positivity of density and pressure for solving com- 
pressible Euler equations using finite volume numerical methods. A general framework is 
established to obtain positivity preserving first and higher order schemes for one and two 
space dimensions. 
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6 Appendix: Positivity of the Lax-Friedrichs Scheme 


We consider in this Appendix, the Lax-Friedrichs scheme for the two dimensional situation 
in Section 3. Fixing a unit vector v, the velocity vector can be written as u = (u,u>) in the 
frame (//, v 1 -). Then, we define 

U = (p, pv, pw, E), E = -p(v 2 + w 2 ) + pe, (6-1) 

H(U) = F(U) ■ v = (pv, pv 2 + p, pvw, ( E + p)v) (6.2) 

The Lax-Friedrichs flux is 

h(V, U) = l -(H(U) + H(V) - fi(V - U)) (6.3) 

with 

P = ll(l U l + C )lloo ( 6 - 4 ^ 

Let us consider UJ, f/" ±1 with positive density and pressure, and set 

up' = (/; - a - M c”, (/;.,)] (6.5) 

Then also has positive density and pressure. An easy way to see that (Sanders [8]) 

consists of introducing a splitting of the equation 

Ut + H(U) X = 0 (6.6) 


by 


U t + (H{U) + (3U) X = 0 (6.7) 

U t + {H(U)~ (W) x = 0 ( 6 . 8 ) 


In other words, we write an exact solver for (5.7)-(5.8), which is easy because there are no 
longer wave interactions with the choice of /3 in (6.4): 

U = Uj- A [//((/") + (HI } 1 - //(!/;_, ) - ( 6 - 9 ) 

b = Uj - \ [//(f/jv, ) - - //(t/;) + /^ 7 ;] ( 6 - 10 ) 
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Now U? + } = is indeed the solution to the Lax-Friedrichs scheme (6.5), and since the 
exact solver preserves positivity (notice that the addition of ±0U is unessential for positivity 
by Galilean invariance), by concave combination (7” +1 satisfies also the positivity property. 
We need the CFL condition 

a ||(M + <011^ < - (6-ii) 

in order to avoid the interaction of waves when solving (5.7)-(5.8) by (5.9)-(5.10). 

Another version of Lax-Friedrichs scheme consists of introducing exact solvers on a stag- 
gered grid. It will also obviously satisfy the positivity property. Both of these Lax-Friedrichs 
schemes in addition satisfy the maximum principle on the specific entropy (Tadmor [10], [5]). 

We would like to conclude this Appendix by a remark raised in [7] on two dimensional 
schemes. Usually two dimensional schemes have to be written under the form (3.4) where the 
approximate solver h(U, V, v ) is indeed a function of three parameters U, V and //, satisfying 
h(U, U, is) = F(U) ■ v. In [7], the notion of “genuinely multidimensional solver” is introduced, 
where the approximate solver is indeed under the special form h(U , V)-u satisfying h{U, U) — 
F(U). The Lax-Friedrichs scheme is obviously not “genuinely multidimensional” because its 
value really depends on the frame (is, v L ) used to write the four components system (6.1)- 
( 6 . 2 ). 

The only example we know of a “genuinely multidimensional” solver is a particular class 
of Boltzmann solvers introduced in [7]. 
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